A Monte-Carlo study of meanders 
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We study the statistics of meanders, i.e. configurations of 
a road crossing a river through n bridges, and possibly wind- 
ing around the source, as a toy model for compact folding of 
polymers. We introduce a Monte-Carlo method which allows 
us to simulate large meanders up to n = 400. By perform- 
ing large n extrapolations, we give asymptotic estimates of 
the connectivity per bridge R — 3.5018(3), the configuration 
exponent 7 — 2.056(10), the winding exponent v = 0.518(2) 
and other quantities describing the shape of meanders. 
Keywords : folding, meanders, Monte-Carlo, tree 
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I. INTRODUCTION 

The concept of folding has an important place in poly- 
mer physics KM- Considering a polymer chain made 
of n identical constituents (the monomers) , the entropy 
of such a system can be obtained by counting the num- 
ber of inequivalent ways of folding the chain onto itself. 
If the model of polymer does not take self-avoidance 
into account, it is then equivalent to the well-known 
Brownian motion. Several more involved models have 
been proposed which study the compact folding of a self- 
avoiding polymer as a Hamiltonian cycle (i.e. a closed, 
self- avoiding walk which visits each vertex) on a regular 
lattice y,Ql. They can also be defined for some kinds 
of random lattices pM] , where each configuration is now 
described by a system of non-intersecting arches which 
connect the pairs of monomers, which are neighbors in 
the real space. 

In the present paper, the compact folding of a polymer 
chain is modeled by a folded strip of stamps, with a com- 
plete pihng of the strip on top of one stamp |^ . It is then 
equivalent to another model of non-intersecting arches, 
the so-called meander problem, which can be summa- 
rized by this simple question: in how many ways M„ can 
a road cross a river through n bridges, and possibly wind 
around the source. 

A related problem can be defined by forbidding the 
winding around the source (i.e. the river is infinite at 
the both ends). It is now equivalent to enumerate the 
"simple alternating transit mazes" M of depth n; it was 
also investigated in connection with Hilbert 16'th prob- 
lem, namely the enumeration of ovals of planar algebraic 
curves M- 



By analogy with some models of statistical mechanics 
like random walks or self-avoiding walks, the meanders 
can be described in the language of critical phenomena. 
In particular, the asymptotic behavior of meanders when 
n is large can be characterized by "critical" exponents. 
However the exact enumeration of meanders is particu- 
larly complicated : there is no known formula for M„ in 
terms of n. By generating all possible configurations, by 
hand or with a computer, the beginning of the sequence 
M„ can be computed |10-15| exactly. As M„ increases 
exponentially with n, the limits of computers are reached 
for n '^ 30 and the estimates of the exponents are too in- 
accurate to validate (or invalidate) some conjectures. 

One should mention that several exact results, for ar- 
bitrarily large n, which are unfortunately not helpful 
to determine the values of the exponents, have been 
obtained with other techniques: random matrix model 
methods ||l^,|l6|-|8[ and an algebraic approach using the 
Temperley-Lieb algebra P3,p0| or the Hecke algebra pi| . 

Many models in statistical physics can be studied by 
Monte-Carlo (or stochastic) methods. With these algo- 
rithms, only a small set of configurations among all the 
possible ones are generated. In principle, the expectation 
of physical quantities (like energy or magnetization) with 
a given law of probability (like the Boltzmann law involv- 
ing an external temperature) can be estimated from these 
randomly generated samples, if their probabilities of gen- 
eration are known. For example, with the Metropolis 
algorithm for classical spin systems [g^li the probability 
of generation is built to be equal to the Boltzmann law 
and the average is done over the generated configurations 
with equal weights. To bypass some difficulties (for ex- 
ample when the phase space has many local minima with 
high free energy barriers between them), it is possible, in 
principal, to generate the random configurations with an- 
other more adapted law and to correct this bias when the 
average is done p^ . 

But for the meanders problem, the situation is quite 
different : the phase space is not easy to built because 
the number A/„ of configurations is unknown and the 
only known efficient method to draw a meander of size 
n is a recurrence over n. Moreover the naive way to use 
this recurrence gives a distribution of meanders which 
is not flat, and this default increases exponentially with 
n. This paper presents a Monte-Carlo method which ex- 
plores the meanders with an almost flat distribution law. 
Furthermore the bias is known and can be corrected ex- 
actly. Therefore, the average can be done over meanders 





FIG. 1. The M4 = 4 meanders of size 4. The road is the 
non-self-intersecting loop. The semi-infinite river is the solid 
half- line, starting at the source (black dot). The size is the 
number of bridges. The winding number w is the number of 
arches crossing the dashed half-line on the right side. The 
upper meanders have no winding {w — 0), but the lower have 
w = 2. 



with equal probabilities. In particular, better estimates 
of critical exponents are obtained. 

After Section 2 devoted to the definitions, and Sec- 
tion 3 which explains the building of meanders by recur- 
rence. Section 4 of this paper describes the Monte-Carlo 
method. The results are presented and discussed in Sec- 
tion 5. 



II. DEFINITIONS 

A meander of size n is a planar configuration of a non- 
self-intersecting loop (road) crossing a half line (semi- 
infinite river with a source) through n points (bridges). 
Two meanders are considered as equivalent if their roads 
can be continuously deformed into each other, keeping 
the bridges fixed: this is therefore a topological equiva- 
lence. We call an arch each section of road between two 
consecutive bridges. So a meander of size n has n bridges 
and n arches. 

The number of different meanders of size n is denoted 
by M„. For example. Mi = 1, M2 = 1, M3 = 2, M4 = 4. 
In Fig. H, the 4 meanders of size 4 are drawn. The M„'s, 
up to n = 29, can be found in Ref. |14| , |l5| . 



In previous articles [|13 14 19f|, these objects were called 



semi- meanders, to distinguish them from the case where 
the line is infinite (river without source). In this paper, 
the river is always a half-line and the word meanders is 
used for convenience. 

We can define ||l^,|lj,|l^ meanders with k connected 
components, i.e. made of one river and k non-intersecting 
roads. But, in this work, we do not include this gener- 
alization and we keep fc = 1. However the Monte-Carlo 
method, used in this article, can be adapted without dif- 
ficulties to an arbitrary fixed fc, and even for varying fc 
with a fugacity q'^. 

As explained with many details in Ref. |lj] , the mean- 
der problem is absolutely equivalent to the problem of the 
compact folding of a strip of stamps because each mean- 
der of size n can be continuously deformed in such a way 



that the "road" becomes a vertical line and the "river" 
becomes a folded strip of n — 1 stamps. We prefer to 
present our results with the meander representation be- 
cause the main recursion relation, described later, seems 
more "natural" in this picture. 

The meander problem has certain similarities with two- 
dimensional self-avoiding walks: a meander is obtained 
by intersecting a closed self-avoiding walk by a half-line 
and keeping only the topological aspect. By analogy, it 
is expected nm that 
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where the estimates given in Ref. |lj| are R = 3.50(1) 
and 7 ~ 2. 

The connectivity R can be reinterpreted as the average 
number of ways of adding a bridge close to the source by 
deforming an arch of a given large meander. Then ln(i?) 
is the entropy per bridge. The configuration exponent 
7 is sensitive to the boundary conditions, for example 
whether the road is closed or open, whether the river is 
infinite or semi-infinite, straight or forked, whether the 
meander is drawn on a planar surface, a sphere or a sur- 
face with higher genus. Conversely, we expect that R 
remains the same for all these boundary conditions. 

It is similar in on-lattice self-avoiding walks problem 
where the connectivity depends on the type of lattice 
(square, honeycomb . . . ) and not on the boundary con- 
ditions, whereas the "universal" configuration exponent 
depends on the boundary conditions, but is not sensi- 
tive to the small scale details of the lattice. For these 
reasons, we think that the numerical value of R is valid 
only for this particular model of meanders. But 7 is ex- 
pected to be more "universal" and to keep its value in 
other variants of the meander problem. Unfortunately, 
the numerical determination of 7 is less precise than R, 
because n"^ describes the correction to the leading expo- 
nential asymptotic behavior i?". 

For a given meander to, the winding number w{m) of 
the road around the source of the river can be defined 
as the minimal number of intersections between the road 
and a half (semi-infinite) line starting at the source and 
extending the river on the opposite side. For an example, 
see Fig. |l[ We define Wn as the average of the winding 
number 



7n—l 



(2) 



over all the meanders to of size n. We can see the winding 
number as the topological end-to-end distance between 
the source (right end of the river) and the infinite (left 
end of the river) . Here the distance between two points 
is simply the minimal number of arches which must be 
crossed to go from one point to the other. By analogy 
with the end-to-end exponent of self-avoiding walk, we 
expect that 




FIG. 2. The height h{i) (resp. h{—i)) is the number of 
arches over (resp. below) the segment i. For this meander 
of size 5, {h{i)} = {0, 1, 2, 3, 2, 1, 0, 1, 2, 1, 0} for i = -5, . . . 5. 
The area is ^ = 13. 

Wn ^ n , (3) 

where the estnuate given in Ref. ||lj] is v = 0.52(1). 

If we study the meanders by leaving free the number k 
of connected components, the problem is equivalent |lj] 
to a random walk on a semi-infinite line and can be stud- 
ied with usual methods of combinatorics. In particular, 
7 = 3/2, R = A and v — 1/2 is the Brownian exponent. 
But, by fixing fc = 1, the problem is drastically more dif- 
ficult and the above values are, to our knowledge, not yet 
known exactly. 

For a given meander m of size n, we label by j = 0, . . . n 
each segment of river in-between two consecutive bridges, 
from right to left. So the rightmost segment (with source) 
is labeled 0, and the leftmost (semi-infinite) segment is 
labeled n. We define the height h(i, m) as the number of 
arches over the segment i, and h{~i,m) as the number 
of arches helow the segment i. An example is given in 
Fig. |. 

For the case i = 0, the both definitions h{+Q,m) 
and h{—Q,m) are equivalent and equal to the wind- 
ing number: h{0,m) = w{m). From the definition, 
we have h{n,m) = h{—n,m) = 0, h{i,m) > and 
h{i + l,m) = h{i, m) ± 1. 

For a given meander m of size n, we define the area 
A{m) as 



A[m) = 2, h{i,m). 



(4) 



For meanders of size n, it can be proved that the max- 
imal area is (n — 1)^ -I- 1 for the two meanders with a 
snail shape (where {^(i)} — 0, 1, 2, .... n — 2, n — 1, n — 
2, . . . , 2, 1, 0, 1, 0} plus the symmetric meander), and the 
minimal area is 2n — 2 (resp. 2n— 1) when n is even (resp. 
odd) for the 2"/^^^ (resp. 2^"^^^/^) snake shaped mean- 
ders characterized by h{i) -I- h{—i) = 2 for < i < n. As 
in the case of the winding number, we will consider the 
average profile height 



hn{i) 



and the average area 
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(a) 





(b) 




FIG. 3. A meander of size n-|- 1 is built from a meander of 
size n with one labeled exterior arch by the following process, 
(a) Add a new bridge on the left side of the river. Cut the 
labeled arch. Stretch its two free ends, (b) Close the arch 
on the opposite side by crossing the new bridge (possibly by 
bypassing the source on the right). The inverse process is the 
following, (b) Open the road at the place of the left most 
bridge, (a) Pull the two free ends and close them on the 
opposite side to form a exterior arch. 



over all the meanders m of size 



III. RECURSION RELATION FOR MEANDERS 

In this section, we describe a recursive algorithm to 
enumerate and built all meanders of a given size n. 
Though it was described in Ref. |0,H, we prefer to 
recall it in details, because the Monte-Carlo method is 
based on this recursion. 

We have different ways to built a meander of size n -I- 1 , 
starting from a meander of size n. Our method consists 
in adding a bridge on the left most part of the river (op- 
posite to the source) and changing the road to cross this 
new bridge. To keep this change minimal, only an exte- 
rior arch is modified (an arch is exterior when no other 
arch surrounds it). 

Take a meander of size n (the parent) and choose (or 
label) one of its exterior arch. By the process described 
in Fig. y, a meander of size n + 1 (the child) is built. A 
parent has as many different children as exterior arches. 
By inverting this process, it appears that each meander 
of size n+1 has one and only one parent. More precisely, 
it is a one-to-one mapping between the meanders of size 
n+1, and the meanders of size n with one labeled exterior 
arch. 

The starting point of the recursion is the unique me- 
ander of size 1. By n — I successive applications of the 



up-down 
symmetry 




FIG. 4. The tree of meanders up to n = 5. For n > 4, only 
a half of the tree is drawn by using the up-down symmetry. 
Then Mi = 1, M2 = 1, M3 = 2, M4 = 4 and M3 = 10. Each 
arrow represents a process as described in Fig. H. 



recursion process, every meanders of size n can be built. 
As shown in Fig. ^ the set of meanders is organized as a 
tree. The root, at level 1, is the starting meander n = 1. 
Each branch between a node on level n and a node on 
level n+1 represents a relation between a parent of size n 
and its child of size n + 1. Apart from n ~ 1, a meander 
(or node) has several exterior arches, then several chil- 
dren (or branches). Their number depends on the precise 
shape of the parent and varies between 2 and n/2 + 1. 

If we want to exactly enumerate the M„ meanders of 
size n, the only method we know is to built and investi- 
gate the tree up to the level n. In particular, we have not 
found a direct recursion between the numbers M„. The 
number of children of each meander has a distribution 
which seems to be erratic and the only way to know it is 
the examination of its shape. Then, to compute M„, the 
work is proportional to Mn. As the M„'s increase expo- 
nentially, the limits of the capabilities of the computers 
are rapidly reached. 

In Ref . H,|lHl , the meanders numbers Af„ up to n = 29 
are given. With a recent computer and more tricks of 
programming, it is perhaps possible to obtain n = 31 
or 32. If the power of computers continues to increase 
exponentially, the best we can expect with the full enu- 
meration without major improvement, is a linear growth 
in n, with a rate of only one new size every 2 or 3 years. 
We have the feeling that this progress is too small to 
change significantly our understanding of the meanders 
problem. 

In this article, we will investigate larger n with a 
Monte-Carlo method increasing more slowly than an ex- 
ponential. But, this stochastic method gives results with 
error bars and the interpretation is more delicate. 



IV. THE MONTE-CARLO METHOD 

As explained in the previous section, the exponential 
growth of the computations with an exact enumeration 
method limits the size of meanders around n ~ 30. To 
reach bigger n, it is then natural to try to study this 
problem with a Monte-Carlo (or stochastic) method. 

As the set of the M„ meanders of a given size n is 
too large to be fully explored, the general idea is to ran- 
domly select a small subset. Then, the measurements are 
done and averaged on the selected meanders. It gives an 
estimator of the exact (but unknown) result, with a un- 
known error. This error has two components : statistical 
fluctuations and bias. 

The statistical fluctuations can be reduced by indepen- 
dently repeating the procedure many times. Then, we 
obtain a histogram of the estimator, with an average and 
a variance. Under the hypothesis of finite variance, the 
statistical fluctuations of the average can be estimated 
by usual formulas of statistics. 

The bias is the difference between the exact result and 
the mathematical expectation of the estimator. If it can 
be exactly calculated, we subtract it from the estimator. 
But, in general, an unknown part remains, which can not 
be reduced by a better statistics. As explained below, by 
adjusting parameters of the simulation, the bias can be 
reduced to become smaller than the statistical fluctua- 
tions. 

In this section, we will first introduce the simplest algo- 
rithm, the one-squirrel method. We will see that its sta- 
tistical fiuctuations grow exponentially with n and they 
are too big for n « 30. Then, we present an algorithm, 
the multi-squirrel method, for which the fluctuations in- 
crease less rapidly. 



A. The one-squirrel method 

The method is based on the recursion relation (see 
Fig. H) , with which the set of meanders is organized as a 
tree (see Fig. ||). The Monte-Carlo squirrel has the fol- 
lowing stochastic behavior. It starts at the root of the 
tree (the meander of size n — 1). It climbs into the tree. 
At each level n, it stands on a node and makes some mea- 
surements concerning the meander of size n, represented 
by this node. Then the squirrel goes to the level n -\- 1 
by choosing at random one of the 6„ branches, starting 
from this node. The squirrel stops at a prefixed level 
n — TT-max- This proccss constitutes one simulation. 

The probability that the squirrel reaches a given me- 
ander of size n is 1/ n;<n ^i- This probability law is not 
flat because the sequence of bi depends on the visited 
nodes. As seen in Fig. 0, for n = 5, the two meanders 
on the left side have a probability 1/8 and the three on 
the right side have 1/12. So, to correct this bias between 
meanders, the squirrel has a weight 



qn 



1=1 



(7) 



which is calculated during its climbing. 

By noting (•) the mathematical expectation (over all 
possible simulations), it is then obvious that (g„) — Mn, 
because the sum runs over M„ possible paths and the 
contribution of a given path is 9„ with probability 1/qn- 

More generally, for some quantity Z (for example the 
winding number), we wish to determine 



E 



Z{m) 



(8) 



where Z{m) is the value of Z on the ?7i-th meander of size 
n. Each simulation gives, at level n, a measurement Z{s) 
on the meander s reached by the squirrel and {qn-Z{s)) = 
Z which is the generalization of {qn) = Mn obtained with 
Z = 1. Then 



Z{s) 



(9) 



is an unbiased estimator of Z (i.e. (z) = Z). 

As usual in Monte-Carlo methods, several simulations 
are made independently and we hope that the average z 
of all the measurements z is close to Z. Unfortunately, 
this method does not work in practice, because the weight 
qn is the product of bi. Although the distribution of each 
bi is regular, the product of many random variables is 
not self-averaging. 

As the sum of ln(6;) is self-averaging (i.e. the observed 
result is closed to its mathematical expectation when n 
goes to infinite) , most of the observed qn are not closed 
to (qn) and 



qn (observed) 



exp^(ln(&,)-(ln6,)) 



(10) 



Kn 



increases like an exponential. Then the averages with 
weight q„ are dominated by exponentially rare events and 
the statistical fluctuations become large. To keep the 
observed average close to the mathematical expectation, 
the number of simulations must increase exponentially 
with n and fluctuations become too big for n '^ 30 or 35. 
As the difficulties increase exponentially with n (as for 
exact enumeration), is is useless to increase the power of 
the computer. We need a new algorithm. 



B. Multi-squirrels method 

We generalize the one-squirrel method, but now with 
a population of S squirrels, which reproduce and die ; 
5 is a fixed parameter during all the simulations. It 
is more simple to choose 5" as S* = Mn„, with at the 
starting point, a squirrel staying at each node of level 
no (meanders of size uq). In this work, rig — 17 and we 



use the up-down symmetry to reduce the population to 
S = Mn/2 = 1664094 squirrels. 

The population evolves from level n to level n -I- 1 by 
the following process. Each squirrel i lives on a node 
Si on level n, connected to bi nodes on level n -I- 1. It 
reproduces and has bi children and each child lives on 
each one of these bi nodes. The total number of chil- 
dren S' = X]i=i ^i is calculated. The ratio i?„ = S' / S is 
an estimate of M„_|_i/M„: the average of the number of 
children per parent. To prevent an exponential growth of 
the population and of the needed computer memory and 
time, the total population is keep constant by decimat- 
ing the children : only S among the S" children survive. 
The choice is made at random with uniform distribution. 
Then the probability of surviving is l/Bn- This decima- 
tion is the single Monte-Carlo step of the algorithm. 

This process is iterated up to reach a prefixed level 
n = Umax'- it gives one simulation. Many independent 
simulations are done and averaged. 

The particular case S = 1 gives the previous method 
with one squirrel. But, as for S' = 1, for every value of S*, 
the probability that a given meander is reached, is not 
uniform. The nodes with small number of "brothers" or 
"cousins" have always a small advantage. But this bias 
becomes smaller when the population S is large. That is 
the main improvement of this method. The limit S* = oo 
corresponds to the exact enumeration. 

To correct the bias, each simulation has a weight 



n^^ 



(11) 



l=no 



and the averages runs over all the simulations with their 
weight. More exactly, for some quantity Z, by keeping 
the notations of Eq. (0) , one simulation with S squirrels 
gives S measurements {Z{si)} for i e [1, S] with a weight 
Qn and the estimator 



qn 



E^(^«) 



is unbiased, i.e. 



(z) = Z. 



(12) 



(13) 



We note that the case Z = 1 gives S ■ {qn) = Mn- 

In order to prove Eq. (p^), we define the operator Sm 
characterizing a given meander m by (5,„(?n') = 1 when 
m = m' and otherwise. Then every operator Z can 
be split up into Z = ^^ Z{m) ■ Sm- As the expectation 
of a sum is always the sum of expectations, we have to 
prove Eq. (^_3|) for the operators Sm only, which becomes 
{qn ■ Am) — 1, where A„i = 1 if the m-th meander is 
occupied by a squirrel and otherwise. Let p represent 
the parent of m in the tree at level n—l. The probability 
that A„j = 1 (i.e. m is occupied) is the product of the 
probability Ap that p was occupied at the level n—l and 



the probability l/i?„_i that its child m survives after the 
decimation process (n — 1 —^ n). 

By using Eg. (pi]) and averaging on the random deci- 
mation (n — 1 — > n) only, 



(A™ 11 Bi) - (Ap J] Bi 



(14) 



l=no 



l=no 



It is a recursion relation between a given meander of size 
n and its parent of size n— 1. By iterating, we go down to 
the ancestor at the starting level ng for which A = 1 and 
the empty product of i?; is 1. It proves that ((?„ • A,„) = 1 
for every meander m of size n, and Eq. ( |l3|) is valid for 
every operator Z. 

As usual, many simulations are done and the measure- 
ments are averaged with their respective weight. A priori, 
it seems that this method has the same defect as the one- 
squirrel method because the weight g„ is the product of 
many Bi , not self-averaging when n becomes large. The 
ratio 
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exp 
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{ln{Bi)-{\nBi)) 



(15) 



between the mathematical expectation and the most fre- 
quently observed (/„ increases like an exponential. But, 
the main improvement of the 22iu7ti-squirrel method is 
that the distribution of Bi becomes narrow when the 
population S of squirrels is large. As Bi = S' /S, the 
fluctuations of Bi are of order 0{l/^/S) because the 
number S" — J^i ^i of children is the sum of S ran- 
dom variables. A Taylor expansion of InBi shows that 
\n{Bi)-{\nBi)=0{l/S) 

Then, with these simple arguments, we can hope that 
the ratio (|l^) grows like 1 -f 0{n/S) and that problems 
appear only when n becomes on the same order than S. 
In our simulations, we observe that the fluctuations grow 
with n faster than this optimistic prediction 0{n/S). In 
fact, the Bis are not independent and the exponential 
function accentuates all deviations. Then we supervised 
carefully the distribution of g„. When n is small, we 
see a regular bell-shaped curve. But, when n increases, 
the distribution becomes asymmetric, with a long and 
irregular tail for the large qn- 

For example, for n — 400, with S = 1664094 squirrels, 
the width a of the distribution is only 12 %, but we ob- 
served rare events with a value of g„ as big as three times 
the average. However, in this case, their contribution to 
the average and fluctuations is not yet problematic. But, 
if we let n increase without control, rare events will dom- 
inate and the results will become hazardous. 

How to choose n and S ? The naive point of view is 
to take n as bigger as possible. But, to make Ng inde- 
pendent simulations with 5* squirrels of size n, the need 
for computer memory is of order 0{n ■ S) and the need 
for computer time is of order 0{n^.S.Ns). If S is large 
enough, the fluctuations are Gaussian and the error bars 



are of order 0{l/^/S.Ns). As n is always limited, we will 
extrapolate to study the asymptotic behavior. For that, 
it is of no help to have large values of n if the error bars 
are too big. So for a fixed computer time, we prefer ac- 
cumulate good statistics by limiting n < 400. Finally for 
a fixed product S.Ng, we prefer to take S — 1664094 as 
bigger as permitted by the memory computer to avoid 
the problem of rare but large fluctuations. 



C. Bias for non-linear observables 

In the previous section, we have seen how to obtain 
unbiased Monte-Carlo estimates of the sum Z over all 
the meanders of size n of some quantity Z (see Eq. (g) 
and its notations). However we are more interested by 
the average Z/M„ over all the meanders of size n. For 
example, the average winding number w„ (see Eq. H) is 
obtained when Z counts the winding. To evaluate R of 
Eq. (1^), we can analyze Af„+i/M„; in this case, Z counts 
the exterior arches. More generally, we want to use non- 
linear combinations of Z and M„ . 

With our Monte-Carlo method, we have seen that one 
simulation gives a measurement z which is an unbiased 
estimator: (z) = Z. With N^ independent simulations, 
we call z the usual average of the Ns measurements z; 
its fluctuations are vNg times smaller. The bar over the 
symbols distinguishes the average of observed values by 
Monte-Carlo method, from the (unknown) mathematical 
expectation, marked with (...). The same work can be 
done with g„ which is a unbiased estimator of M„. 

We must be careful with the Monte-Carlo estimate of 
Z/Mn- For example, the average of the ratio z/qn gives 
bad results. It is better to compute the ratio of the av- 
erages z/qn- Indeed, with a Taylor expansion of z and 
qn around their mathematical expectations Z and M„, 
the bias (defined as the difference of the mathematical 
expectation {z/qn) and the target Z/Mn) 



{z/qn)-Z/Mn^O{l/Ns). 



(16) 



It can be neglected if it is smaller than the stochastic 
fluctuations. Usually, in Monte-Carlo simulations, this 
problem disappears because several millions of indepen- 
dent measurements are done. But, in this work, the sit- 
uation is quite different. In fact, each simulation is a 
complex process involving millions of squirrels, and the 
number Ns of simulations is small. 

Of course, we can not compute this bias exactly, oth- 
erwise we would have already subtract it from measure- 
ments. But we can estimate it by the following process. 
The set of simulations is divided into Ng/'iP subsets, with 
2^ simulations each. In each subset, z/qn is computed. 
We obtained Ns/2'p independent values, one for each sub- 
set; by usual formulae of statistics, we compute their av- 
erage E'^p'^ and the error bars. This work is done for all 
integer p between p — Q (each subset contains only one 



log2 N, 



(only one set with all the Ng 



simulation) and p 
simulations) . 

Which value of p is the best ? Following Eq. (16), 
the bias of E'^p'' is expected to decrease like 1/2^. For 
small values of p, we observe really a dependency of i^^^-* 
on p: the bias is visible. But, for p > 5, variations be- 
come smaller than statistical error bars: the size 2^ of 
subsets is sufficiently large to neglect the bias. But, if 
p is close to its maximum, the number of subsets be- 
comes very small and the error bars are not properly es- 
timated. As in our work Ng — 8192 = 2-'^'^, we finally keep 
p — 7: the estimator E^'^' is computed with 64 indepen- 
dent subsets of 128 simulations each. This method was 
used for all quantities presented below. It is valid, not 
only for ratios like Z/Mn, but also for non- linear func- 
tions like ln(M„4.i/M„). It could also be possible to use 
more complex estimators. For example, the combination 
2^(p) _ £'(p-i) cancels the order 1/2^ of the bias. 



V. RESULTS 

In this section, we describe the results obtained by 
our Monte-Carlo multi-squirrels method. After several 
tests, we used a population of 5 = 1664094(= A/17/2) 
squirrels for meanders with size up to n = 400. To do 
Ng = 8192 independent simulations, we have used during 
8 days a parallel computer (the Cray T3E of the Cea- 
Grenoble) with 128 processors (Dec- alpha at 375 MHz) 
and 13 Gigabytes of total memory, equivalent to 24000 
hours of workstation cpu time. 

We have verified that the results are stable when 5' 
(the population) increases. More exactly, the tests with 
smaller S have larger error bars, but are compatible with 
results and error bars obtain ed with the largest S. As 
explained in the Sect. 



[VB, we have carefully checked 



that S = 1664094 is sufficiently large to explore sizes of 
meander up to n = 400. 



A. Enumeration 

We want to measure R and 7, which describe the 
asymptotic behavior of the number of meanders M„ ~ 
c R'"'/n^ for large n. The entropy Ini? can be estimated 
by ln(M„/M„_i). But it appears that the sub-sequences 
M2n and M2n+i have an alternating sub- leading correc- 
tion. We have estimated it to be u(— l)"/(n Inn) with 
u = 0.5(1). This alternating effect is dramatically ampli- 
fied by the ratio Mn/Mn-i- So it is better to consider 



Lr. 



1 



■ln(Af„/M„_2), 



(17) 



with a jump from n ~ 2 to n. But even with this pre- 
caution, the reader can still see on the following figures 
a small parity effect. To estimate Ln, we have used the 
procedure described in Sect. IV C. 
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FIG. 5. ln_R and 7 : plot of the Monte-Carlo estimate of 
Ln — In 3.5 -|- 2/n, for n between 50 and 400, versus 1/n. The 
limit when x goes to is In(_R/3.5), and the (negative) slope 
is 2 — 7. The error bars are not drawn ; their maximum is 
10~^, then they are smaller than the symbols. A parity effect, 
between the odd and even n, is visible. A linear extrapolation 
gives R = 3.5019(2) and 7 = 2.056(10). 
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FIG. 6. InR and a under the hypothesis 7 = 2: the 
same plot as FigJS but the x-axis is l/(n In n). The (negative) 
slope is —a. A linear extrapolation gives R = 3.5017(2) and 
a = 0.25(5). 



As we expect L„ ~ Ini? — 7/n for large n, by plot- 
ting y = Ln versus x = 1/n, we can estimate Ini? 
(limit when x goes to 0) and the exponent 7 (asymp- 
totic slope). In Fig. ^ we have plotted the Monte-Carlo 
estimate of L„ — In 3.5 -I- 2/n versus 1/n for n between 
50 and 400. We have arbitrarily subtracted the linear 
function y = In 3.5 — 2x, to reduce the amplitude of y; 
we obtain a figure where the small quantities 2 — 7 (re- 
maining slope), ln(i?/3.5) (limit when x goes to 0) and 
curvature (deviation to the expected linear behavior) are 
more visible. The curvature remains small and a linear 
extrapolation gives a limit between 0.0005 and 0.0006 
with an estimated slope 0.056(10). Then 



i? = 3.5019(2) and 7 = 2.056(10). 



(18) 



With the assumption that the asymptotic behavior is 
Mn ^ c R^ /rf , the conjecture 7 = 2 is incompatible 
with these simulations. But, we can try another asymp- 
totic shape, for example 



M^r^C—r 



1 



rf In" n ' 



(19) 



by introducing a new exponent a. In Fig. o, we have 
plotted Ln — In 3.5 -I- 2/n (as in Fig. ^), but now versus 
l/(nlnn). With this transformation of the x-axis, a lin- 
ear behavior corresponds to 7 = 2 and the slope measures 
—a. A linear extrapolation gives 



i? = 3.5017(2), 7 = 2 and 



0.25(5). (20) 



with a compatible with the simple fraction 1/4. 

How shall we choose between both results Eq. ( p^ or 
Eq. (20) ? We notice that the quality of the alignment of 
points is the same in Fig. || and Fig. |^. In fact, it is very 
difficult to distinguish between a logarithmic law and a 
power law with such a small exponent, with statistical 
error bars and when the amplitude of In n is small. In 
fact, for n close to n', the term Q;/(n Inn) looks like a'/n 
with a' / a ~ \/\nn' -I- 1/ln n'. Then, for any choice 
of a not too large, 7 = 2.056 — 0.22a gives a class of 
acceptable behaviors. 

We have tried to use more sophisticated extrapolation 
methods. For example, with a fixed jump i, (ni„ — (n — 
i)Ln-i)/i gives theoretically the same limit ln_R but by 
removing the term 7/n. Thus •n?{Ln — Ln-i)/i gives a 
direct estimate of 7. Unfortunately these kinds of deriva- 
tive amplify statistical errors, and the results are compat- 
ible but less precise than the previous estimates. 

So, with our numerical simulations, we cannot say if 
a logarithmic factor is present or not. However, R = 
3.5018(3), and we can exclude R = 3.5 and the conjecture 
7 = 2 without logarithmic factor {a — 0). 



B. Winding 

We will present the Monte-Carlo results of the expo- 
nent I', which describe the asymptotic behavior of the av- 
erage winding number u;„ ~ n"^. To avoid problem with 



winding exponent v 
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FIG. 7. Winding exponent f : plot of the Monte-Carlo 
results of ln(w;„ -|- 1) — 1/2 Inn, for n between 8 and 400, 
versus In n. The slope is 0.018; it is a measurement of i/ — 1/2. 
The error bars are not drawn; their maximum is S-IO"**, then 
they are smaller than the symbols. 
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bias, we have used the procedure described in Sect. IV C . 
By plotting ln(ii;„) versus ln(n), the asymptotic slope 
will be a measurement of i^. In Fig. [7|, we have plot- 
ted y = ln(w„ -|- 1) — 1/2 Inn versus x = Inn. We 
have arbitrarily considered ln(w„ -I- 1) and not ln(w„) 
because w„ -I- 1 is less sensitive than Wn to the finite size 
effects |1J]. As the main question is to know whether 
ly = 1/2 or not, we have arbitrarily subtracted the linear 
function y = x/2. Then the variation of y is reduced and 
we obtain a figure where i^ — 1/2 (residual slope) and the 
curvature are more visible. We see that the curvature is 
small and a linear extrapolation gives v = 0.518. As it 
is difficult to estimate the errors with the data of Fig. M, 
we have also tried more sophisticated quantities like 



G,(n) = ^ln(^^;^^±^ 



(21) 



which are discrete derivatives of ln(w„ -I- 1) with step i. 
They give a direct value for v, but unfortunately the sta- 
tistical fluctuations are amplified by this differentiation 
and the uncertainty over ly is of order 0.002. 

We have seen that, for the exponent 7, a behavior 
with logarithmic correction is not excluded by the Monte- 
Carlo data. So we have tried to fit the winding number 
with 



ii/2 in° 



(22) 



With this hypothesis, aplot of y = ln(?i;„+l)— l/21nn, as 
in Fig. n, but now versus a; = In In n would give a straight 
line with slope a. But the curvature is much stronger 
than that of Fig. m So we dismiss this hypothesis and 
conclude that 
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(w+1)/(w„+1) 

FIG. 8. Plot of («;„ + l)P„{'w)/2 versus (w + l)/{w„ + 1) 
for n =30, 40, 50 ... 400. The points accumulate on a scaling 
function /. The error bars are not drawn; their maximum is 
6.10-*. 



ly = 0.518(2) 



(23) 



With the assumption that the asymptotic behavior is 
a simple power-law, the Brownian value v = 1/2 is in- 
compatible with these simulations. 



C. Probability distribution of winding number 

We define the probabihty distribution P„ (w) of wind- 
ing number as the fraction of the meanders of size n with 
w windings. We expect [H the asymptotic scaling be- 
havior 



PnH 



W„ 



1 



/ 



1 



1 



(24) 



with a scaling function f{x), where Wn is the average 
winding number. With the factor 2, the integral of / is 
normalized to 1 because w and n are integers with the 
same parity. In Fig. g, we plot y = (w„ -I- l)P„(u;)/2 
versus x ^ {w + l)/(w„ + 1) for different values of n. 
To define the scaling variable x, we prefer to take w + 1 
instead of w to reduce finite size effects. 

We see that the points accumulate on a smooth curve, 
which represents the scaling function f{x). By analogy 
with the end-to-end distribution for polymers, we ex- 
pect Q a power law behavior, f{x) ^ x^ , for small x, 
and a behavior f{x) ~ exp(— const. x'^), for large x. Our 
data give the estimates 



= 1.7(1) and (5 = 2.3(1) 



(25) 



To obtain a better precision, it would require to have 
larger values of n. 
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FIG. 9. Plot of the scaled height n{hn{i) + 1)/(A„ + 2n) 
versus the scaled coordinate x = i/n. Only two set of data 
(n = 40 and 400) are shown in order not to overload the figure. 
The points accumulate on a scaling function p{x). We can see 
small deviations for a; = and x = 1, which are analyzed in 
the text. 



D. Height and area 

We are interested by the asymptotic behavior, for large 
size n, of the average area A„ (see Eq. ^) and average 
height hn{i) (see Eq. ^). 

The label i is the "horizontal" coordinate and varies 
between —n and n. So, we introduce the scaled variable 
X = i/n, with — 1 < x < 1. We expect [|lj| that 



hninx) 



Ar. 



p{x) 



(26) 



for large n by fixing x, where p{x) is a scaling function 
with integral normalized to 1. In Fig. y, we have plotted 
y — n {hn{i) + l) / {An+2n) versus x = i/n for various size 
n. Only the positive i are shown because hn(i) is symmet- 
ric after summing over all the meanders. More precisely, 
we have plotted the average of {/i„(z) -t- /i„(— «)}/2 over 
the Monte-Carlo samples of meanders, which is equiva- 
lent to the average of hn{i) over these samples, plus those 
obtained by the left-right symmetry (i — > — i). As in the 
previous figures, we prefer to take (/i„ -f- 1) and (A„ + 2ri) 
to reduce finite size effects. 

We see that, for < a; < 1, the points accumulate 
on a smooth curve, which represents the scaling function 
p{x). This shape is not a half-circle, as it would be for a 
random- walk on a semi-infinite line p4| . 

As hn(0) is the winding number Wn for the particular 
case a; = 0, Eq. ( |2q ) can be valid only if hn{nx) scales as 
n"^ with the same exponent v for all x. Consequently, the 
area An scales as n""^^. In a previous section, we have 
numerically determined v = 0.518(2) by extrapolation 
of Wn — hn(0). The same work with /i„(n/2) and An 
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FIG. 10. Log-log plot of the height hn(n — i) versus the 
distance to the boundary i for various sizes n. The points 
accumulate on a hmiting curve, which can be fitted by i'* 
with (^ = 0.64(2). 



gives V — 0.517 for both, compatible with the previous 
estimates, but less precise because the finite size effects 
are stronger. 

In Fig. 0, small deviations appear between the curves 
for n=40 and 400, at the boundary (x ~ ±1) and in the 
middle of the meander (x ~ 0). 

In order to understand why finite-size effects are im- 
portant near the boundary, Fig.O is a plot of ft,„(n — i) 
versus i, with various n. Clearly, when n is large, curves 
accumulate on a limiting curve 



h{i) = lim hn{n — i). 



(27) 



As this function can be fitted by a straight line on this 
"log-log" plot, we define a new exponent (j) by 



hit) 



,•'/' 



with (j) = 0.64(2) 



(28) 



when i is large. As by construction, hn{n) ~ and ft.„(n— 
1) = 1 for all n, /i(0) = and hiV) — 1 satisfy exactly 

If we set i = ny, by using Eq. (ES) valid when i is fixed, 
we obtain hn{n—i) ^ {nyY . On the other hand, by using 
Eq. ( p6|) which is valid when y is fixed, we obtain now 
n'^p{l — y). As the exponents v and </) differs, these two 
regimes are incompatible and the behavior of hn{n — i) 
depends on the order in which n and i go to infinity. In 
particular, the domain of validity of Eq. (^8|) is reduced 
to the single point x = 1 for Eq. ( pGJ ) . That explains why 
the rightmost dots (i ~ n) in Fig. ^ for the small size 
are not superimposed on the curve obtained for the large 
size. 

Here, the exponent 4> is the surface critical exponent, 
while V is the bulk critical exponent. Near the boundary. 
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FIG. 11. Plot of the hn{i) — hn{0) versus i for small i, 
near the middle of meander, for various size n. By symme- 
try, only the positive i are drawn. As the heights h„{i) are 
shifted by /i„(0) and not scaled, we observe effects smaller 
than the asymptotic behavior n" p{Q). The points accumulate 
on a limiting curve h{i). 



h{i) is small and the condition h{i) > limits appreciably 
the ffuctuations toward the bottom. This effect is so 
strong that the exponent is changed and (j) > v. This 
situation is reminiscent of other critical phenomena [g4| , 
like the self-avoiding walk near a surface. Our result is 
to be contrasted with the case of a random walk on a 
semi-infinite line for which the surface exponent keeps 
its Brownian value 1/2. 

Finite-size effects observable near the middle of mean- 
der (x ^ 0) can be explained with Fig. O, which is a 
plot of hn{i) — hn{Q) versus i for small i with various n. 
Clearly, when n is large, curves accumulate on a limiting 
curve 



h{i)^ lim {K{i)~K{0)). 



(29) 



If we make the hypothesis that the behavior of h{i) is 
compatible with Eq. (Eq) by inverting the limits n large 
and X small, the consequences would be that p{x) has a 
cusp at x = with a infinite derivative p{x) ^ p(0) -I- x'^ 
when X is small, and h{i) ^ i^ . But this power law 
behavior of h(i) is not observed. Then the asymptotic 
behavior of h^ii) with i fixed and n large is given by 
Eq. (P6|), i.e. n''p(0), plus finite corrections of order h(i). 

This cusp is an another boundary effect since the point 
i = is the source of the river. Let us consider three 
consecutive heights {/i(« — 1, ttt.), /i(«, ttt.), /i(«-(-l, m)} for a 
given meander TO. By definition, /i(z-|-1,to) = /i(i,TO)±l. 
Then, for a generic i, the couple of its neighbors can have 
4 respective values \h[i—\\ /i(i+l)} = {/i(i)±l, /i(i)±l}. 
But, for the special case i = 0, the situation /i(— 1) = 
h(X) = h{0) — 1 happens only if a single arch connects 
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the first bridge (near the source) to itself, by drawing a 
httle circle around the source, without visiting the other 
bridges. This is forbidden if we insist in having only one 
connected component. In other word, the neighborhood 
of the source limits the fluctuations of h{0) toward the 
top. In particular, for every meander, h{0) < {h{—l) + 
h{l)}/2. That explains why, on average, h{0) < h{l). 

To understand why /i(l) < h{2), it is mandatory to 
consider more complex forbidden situations for {h{—2), 
h{—l), h{0), h{l), h{2)}. For example, a circle connect- 
ing bridges 1 and 2 with an upper and a lower arch, with 
h{±l) = h{±2) + 1, is forbidden. 

More generally, for a given i , the presence of the source 
forbids the systems of arches connecting the bridges 
{1,2,..., i} with a closed road without visiting the other 
bridges {i + 1, . . .}. Qualitatively, it is a repulsive "force" 
which favors the connection of bridge i with bridges j > i 
and gives a concave shape for small i. 

As the forbidden situations arc more and more complex 
when i becomes large, their statistical effects decrease 
and this repulsive force has a finite range. In the end, 
the summation over all forbidden situations gives finally 
this cusp with a finite amplitude described by h(i). 



From a Monte-Carlo point of view, the proposed al- 
gorithm can be used, in principal, for any combinatorial 
problem described by a tree : the essential ingredient is 
to know, for a given node, the number of branches. How- 
ever algorithms are rare in this kind of problems and bet- 
ter variants or other algorithms could be without doubt 
invented. 

About the meanders, with these estimates, it appears 
that the critical exponent are not simple fractions as 1/2 
or 7/2, as conjectured by previous studies |1^. Of course, 
Monte-Carlo simulations cannot determine the exact val- 
ues, but can confirm or invalidate analytical proposals, 
while waiting for a rigorous solution. 
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VI. CONCLUSION 



In this paper, we have presented a Monte-Carlo 
method to investigate a phase space described by a deter- 
ministic but irregular tree (i.e. the number of branches 
at each node is not fixed). With a naive random climbing 
on the tree (the one-squirrel method), the probability of 
a path depends on the number of branches encounter at 
each node. For the meanders problem, the ratio between 
extremal weights increases exponentially with the size : 
consequently the most part of the computer time is de- 
voted to generate configurations with small weight, and 
only a exponentially small number of configurations with 
high weight contribute efficiently to the average. 

With the multi-squiTTel method, the distribution be- 
comes almost flat : the bias, i.e. the ratio between ex- 
tremal weights increases very slowly and never exceeds 3 
in our simulations. Moreover, this bias is exactly known 
during the simulation, then it can be corrected to aver- 
age over all meanders with a uniform distribution. As 
usual with Monte-Carlo simulations, results suffer from 
statistical fluctuations which decrease, in the best case, 
like the square root of the computer time. 

After a simulation on a parallel computer with 3 years 
of cpu time in single processor units, results with small 
errors bars have been obtained for meanders up to size 
n = 400. Under some hypothesis inspired by the analogy 
with random walks problems, large n extrapolation can 
be done for the enumeration (see Eq. Ilq and Fig. @), for 
the distribution (see Eq. GB and Fig . pFand the average 
of the winding number (see Eq. ^ and Fig. [7|) and the 
shape (see Section VD) of meanders. 
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